function OO_rsr2000 = cf_rsr2000(fc_multiplier, vc_multiplier,...
    cost2_0, cost2_1_all, D_input, D, N,...
    D_enroll_st, ...
    alpha, ...
    Delta_xi_all, cost1_e_all, mu_f_all, sigma_f_all,...
    K, pr_y, N_y, resource_y1, resource_y2, y_medi, x_fc, pr_fc, ...
    rho, crra, sigma_c1, beta_c, beta_f, T1, T2, B1_c, B2_c, B1_f, B2_f, delta, ...
    lr_target, N_f, entry_exitflag_vec, id_ist_old, myopic, Er, cost_weight, use_cost_weight, root_dir)

% indicator for whether RSR 2000 in place
rsr2000_yes = zeros(size(D_input.s1,1), 1); % NN x 1
temp1 = isnan(D_input.reg_year);
temp2 = (D_input.year>D_input.reg_year);
rsr2000_yes(temp1==1) = 0;
rsr2000_yes(temp1==0) = temp2(temp1==0); % NN x 1
rsr2000_yes_mat = repmat(rsr2000_yes, 1, size(D_input.p2,2)); % NN x KM

disp(' ')
disp('====================================================')
disp('No RSR 2000 vs. All RSR 2000')

% loop over different regimes
for ss=1:3
    tic0=tic;
    disp(' ')

    if ss==1 % no RSR 2000 regime
        disp('***************************')
        disp('ss=1, No RSR 2000 regime')

        % FC: for rsr2000_yes=1, lower
        cost2_0_cf = cost2_0;
        cost2_0_cf(rsr2000_yes_mat==1) = (1/(1+fc_multiplier))*cost2_0(rsr2000_yes_mat==1);

        % VC: for rsr2000_yes=1, lower
        cost2_1_all_cf = cost2_1_all;
        cost2_1_all_cf(rsr2000_yes_mat==1) = (1/(1+vc_multiplier))*cost2_1_all(rsr2000_yes_mat==1);

    elseif ss==2 % all RSR 2000 regime
        disp('***************************')
        disp('ss=2, All RSR 2000 regime')

        % FC: for rsr2000_yes=0, increase
        cost2_0_cf = cost2_0;
        cost2_0_cf(rsr2000_yes_mat==0) = (1+fc_multiplier)*cost2_0(rsr2000_yes_mat==0);

        % VC: for rsr2000_yes=0, increase
        cost2_1_all_cf = cost2_1_all;
        cost2_1_all_cf(rsr2000_yes_mat==0) = (1+vc_multiplier)*cost2_1_all(rsr2000_yes_mat==0);

    elseif ss==3 % baseline
        disp('***************************')
        disp('ss=3, Baseline regime')

        cost2_0_cf = cost2_0;
        cost2_1_all_cf = cost2_1_all;

    end

    % initial guess on eqm outcomes: use observed data
    p1_ini = D_input.p1; % NN x 1
    p2_ini = D_input.p2; % NN x KM
    n_ini_unmerge = D_input.n_insurers; % NN x 1: 1 for majors, # fringes for the fringe
    n_ini = zeros(N.st,1); % N_st x 1. bring down to market level
    for mm=1:N.st
        n_ini(mm) = unique(n_ini_unmerge(D_input.mkt==mm & D_input.major==0));
    end

    % solve the counterfactual eqm
    [p1_cf, p2_cf, n_cf, ...
        profit_cf, non_converge_all_cf, solve_market_all_cf, fsolve_all_cf, ...
        EV_y_cf, EV_cf, enrollment_cf, enrollment_m_cf, enrollment_f_cf, enrollment_per_m_cf, enrollment_per_f_cf, non_converge_cf, solve_market_cf,...
        s1_mkt_cf, s2_mkt_cf, rs_share_cf]...
        = solve_eq(alpha, ...
        Delta_xi_all, cost1_e_all, cost2_1_all_cf,  mu_f_all, sigma_f_all,...
        cost2_0_cf, D_input, N, ...
        K, pr_y, N_y, resource_y1, resource_y2, y_medi, x_fc, pr_fc, ...
        rho, crra, sigma_c1, beta_c, beta_f, T1, B1_c, B2_c, B1_f, B2_f, delta, ...
        lr_target, N_f, entry_exitflag_vec, id_ist_old, myopic, Er,...
        p1_ini, p2_ini, n_ini, cost_weight, use_cost_weight);

    % save all outcomes
    OO_rsr2000(ss).p1_cf=p1_cf;
    OO_rsr2000(ss).p2_cf=p2_cf;
    OO_rsr2000(ss).n_cf=n_cf;
    OO_rsr2000(ss).profit_cf=profit_cf;
    OO_rsr2000(ss).non_converge_all_cf = non_converge_all_cf;
    OO_rsr2000(ss).solve_market_all_cf = solve_market_all_cf;
    OO_rsr2000(ss).fsolve_all_cf = fsolve_all_cf;
    OO_rsr2000(ss).EV_y_cf = EV_y_cf;
    OO_rsr2000(ss).EV_cf = EV_cf;
    OO_rsr2000(ss).enrollment_cf = enrollment_cf;
    OO_rsr2000(ss).enrollment_m_cf = enrollment_m_cf;
    OO_rsr2000(ss).enrollment_f_cf = enrollment_f_cf;
    OO_rsr2000(ss).enrollment_per_m_cf = enrollment_per_m_cf;
    OO_rsr2000(ss).enrollment_per_f_cf = enrollment_per_f_cf;
    OO_rsr2000(ss).non_converge_cf=non_converge_cf;
    OO_rsr2000(ss).solve_market_cf=solve_market_cf;
    OO_rsr2000(ss).s1_mkt_cf = s1_mkt_cf;
    OO_rsr2000(ss).s2_mkt_cf = s2_mkt_cf;
    OO_rsr2000(ss).rs_share_cf = rs_share_cf;

    toc0 = toc(tic0);
    disp(' ')
    disp(['Time elapsed for current ss (min) = ', num2str(toc0/(60), 5)]);
end
% save
save('cf_rsr2000.mat')
